!> \file GFS_phys_time_vary.scm.F90
!!  Contains code related to GFS physics suite setup (physics part of time_vary_step)

!>\defgroup mod_GFS_phys_time_vary GFS Physics Time Update
!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor,
!! aerosol, IN&CCN and surface properties updates.
!> @{
   module GFS_phys_time_vary

      use machine, only : kind_phys

      use mersenne_twister, only: random_setseed, random_number

      use ozne_def, only : levozp, oz_coeff, oz_lat, oz_pres, oz_time, ozplin
      use ozinterp, only : read_o3data, setindxoz, ozinterpol

      use h2o_def,   only : levh2o, h2o_coeff, h2o_lat, h2o_pres, h2o_time, h2oplin
      use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol

      use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm, iamin, iamax, jamin, jamax
      use aerinterp,  only : read_aerdata, setindxaer, aerinterpol, read_aerdataf

      use iccn_def,   only : ciplin, ccnin, ci_pres
      use iccninterp, only : read_cidata, setindxci, ciinterpol

      use cires_tauamf_data,   only:  cires_indx_ugwp,  read_tau_amf, tau_amf_interp
      use cires_tauamf_data,   only:  tau_limb,  days_limb, ugwp_taulat

      !--- variables needed for calculating 'sncovr'
      use namelist_soilveg, only: salp_data, snupx
      use set_soilveg_mod, only: set_soilveg

      ! --- needed for Noah MP init
      use noahmp_tables, only: laim_table,saim_table,sla_table,      &
                               bexp_table,smcmax_table,smcwlt_table, &
                               dwsat_table,dksat_table,psisat_table, &
                               isurban_table,isbarren_table,         &
                               isice_table,iswater_table

      implicit none

      private

      public GFS_phys_time_vary_init, GFS_phys_time_vary_timestep_init, GFS_phys_time_vary_timestep_finalize, GFS_phys_time_vary_finalize

      logical :: is_initialized = .false.

      real(kind=kind_phys), parameter :: con_hr        =  3600.0_kind_phys
      real(kind=kind_phys), parameter :: con_99        =    99.0_kind_phys
      real(kind=kind_phys), parameter :: con_100       =   100.0_kind_phys
      real(kind=kind_phys), parameter :: missing_value = 9.99e20_kind_phys
      real(kind=kind_phys), parameter :: drythresh     =   1.e-4_kind_phys
      real(kind=kind_phys), parameter :: zero          =     0.0_kind_phys
      real(kind=kind_phys), parameter :: one           =     1.0_kind_phys

      contains

!> \section arg_table_GFS_phys_time_vary_init Argument Table
!! \htmlinclude GFS_phys_time_vary_init.html
!!
!>\section gen_GFS_phys_time_vary_init GFS_phys_time_vary_init General Algorithm
!! @{
      subroutine GFS_phys_time_vary_init (                                                         &
              me, master, ntoz, h2o_phys, iaerclm, iccn, iflip, im, nx, ny, idate, xlat_d, xlon_d, &
              jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl,fhour,          &
              jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm,            &
              jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap,              &
              do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,                            &
              isot, ivegsrc, nlunit, sncovr, sncovr_ice, lsm, lsm_noahmp, lsm_ruc, min_seaice,     &
              fice, landfrac, vtype, weasd, lsoil, zs, dzs, lsnow_lsm_lbound, lsnow_lsm_ubound,    &
              tvxy, tgxy, tahxy, canicexy, canliqxy, eahxy, cmxy, chxy, fwetxy, sneqvoxy, alboldxy,&
              qsnowxy, wslakexy, albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, albdvis_ice,  &
              albdnir_ice, albivis_ice, albinir_ice, emiss_lnd, emiss_ice, taussxy, waxy, wtxy,    &
              zwtxy, xlaixy, xsaixy, lfmassxy, stmassxy, rtmassxy, woodxy, stblcpxy, fastcpxy,     &
              smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy,   &
              slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, lsm_cold_start, nthrds,    &
              errmsg, errflg)

         implicit none

         ! Interface variables
         integer,              intent(in)    :: me, master, ntoz, iccn, iflip, im, nx, ny
         logical,              intent(in)    :: h2o_phys, iaerclm, lsm_cold_start
         integer,              intent(in)    :: idate(:)
         real(kind_phys),      intent(in)    :: fhour
         real(kind_phys),      intent(in)    :: xlat_d(:), xlon_d(:)

         integer,              intent(inout) :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:)
         real(kind_phys),      intent(inout) :: ddy_o3(:),  ddy_h(:)
         real(kind_phys),      intent(in)    :: ozpl(:,:,:), h2opl(:,:,:)
         integer,              intent(inout) :: jindx1_aer(:), jindx2_aer(:), iindx1_aer(:), iindx2_aer(:)
         real(kind_phys),      intent(inout) :: ddy_aer(:), ddx_aer(:)
         real(kind_phys),      intent(in)    :: aer_nm(:,:,:)
         integer,              intent(inout) :: jindx1_ci(:), jindx2_ci(:), iindx1_ci(:), iindx2_ci(:)
         real(kind_phys),      intent(inout) :: ddy_ci(:), ddx_ci(:)
         integer,              intent(inout) :: imap(:), jmap(:)
         logical,              intent(in)    :: do_ugwp_v1
         real(kind_phys),      intent(inout) :: ddy_j1tau(:), ddy_j2tau(:)
         integer,              intent(inout) :: jindx1_tau(:), jindx2_tau(:)

         integer,              intent(in)    :: isot, ivegsrc, nlunit
         real(kind_phys),      intent(inout) :: sncovr(:), sncovr_ice(:)
         integer,              intent(in)    :: lsm, lsm_noahmp, lsm_ruc, vtype(:)
         real(kind_phys),      intent(in)    :: min_seaice, fice(:)
         real(kind_phys),      intent(in)    :: landfrac(:)
         real(kind_phys),      intent(inout) :: weasd(:)

         ! NoahMP - only allocated when NoahMP is used
         integer, intent(in) :: lsoil, lsnow_lsm_lbound, lsnow_lsm_ubound
         real(kind_phys),      intent(in)    :: zs(:)
         real(kind_phys),      intent(in)    :: dzs(:)
         real(kind_phys),      intent(inout) :: tvxy(:)
         real(kind_phys),      intent(inout) :: tgxy(:)
         real(kind_phys),      intent(inout) :: tahxy(:)
         real(kind_phys),      intent(inout) :: canicexy(:)
         real(kind_phys),      intent(inout) :: canliqxy(:)
         real(kind_phys),      intent(inout) :: eahxy(:)
         real(kind_phys),      intent(inout) :: cmxy(:)
         real(kind_phys),      intent(inout) :: chxy(:)
         real(kind_phys),      intent(inout) :: fwetxy(:)
         real(kind_phys),      intent(inout) :: sneqvoxy(:)
         real(kind_phys),      intent(inout) :: alboldxy(:)
         real(kind_phys),      intent(inout) :: qsnowxy(:)
         real(kind_phys),      intent(inout) :: wslakexy(:)
         real(kind_phys),      intent(inout) :: albdvis_lnd(:)
         real(kind_phys),      intent(inout) :: albdnir_lnd(:)
         real(kind_phys),      intent(inout) :: albivis_lnd(:)
         real(kind_phys),      intent(inout) :: albinir_lnd(:)
         real(kind_phys),      intent(inout) :: albdvis_ice(:)
         real(kind_phys),      intent(inout) :: albdnir_ice(:)
         real(kind_phys),      intent(inout) :: albivis_ice(:)
         real(kind_phys),      intent(inout) :: albinir_ice(:)
         real(kind_phys),      intent(inout) :: emiss_lnd(:)
         real(kind_phys),      intent(inout) :: emiss_ice(:)
         real(kind_phys),      intent(inout) :: taussxy(:)
         real(kind_phys),      intent(inout) :: waxy(:)
         real(kind_phys),      intent(inout) :: wtxy(:)
         real(kind_phys),      intent(inout) :: zwtxy(:)
         real(kind_phys),      intent(inout) :: xlaixy(:)
         real(kind_phys),      intent(inout) :: xsaixy(:)
         real(kind_phys),      intent(inout) :: lfmassxy(:)
         real(kind_phys),      intent(inout) :: stmassxy(:)
         real(kind_phys),      intent(inout) :: rtmassxy(:)
         real(kind_phys),      intent(inout) :: woodxy(:)
         real(kind_phys),      intent(inout) :: stblcpxy(:)
         real(kind_phys),      intent(inout) :: fastcpxy(:)
         real(kind_phys),      intent(inout) :: smcwtdxy(:)
         real(kind_phys),      intent(inout) :: deeprechxy(:)
         real(kind_phys),      intent(inout) :: rechxy(:)
         real(kind_phys),      intent(inout) :: snowxy(:)
         real(kind_phys),      intent(inout) :: snicexy(:,lsnow_lsm_lbound:)
         real(kind_phys),      intent(inout) :: snliqxy(:,lsnow_lsm_lbound:)
         real(kind_phys),      intent(inout) :: tsnoxy (:,lsnow_lsm_lbound:)
         real(kind_phys),      intent(inout) :: smoiseq(:,:)
         real(kind_phys),      intent(inout) :: zsnsoxy(:,lsnow_lsm_lbound:)
         real(kind_phys),      intent(inout) :: slc(:,:)
         real(kind_phys),      intent(inout) :: smc(:,:)
         real(kind_phys),      intent(inout) :: stc(:,:)
         real(kind_phys),      intent(in)    :: tsfcl(:)
         real(kind_phys),      intent(in)    :: snowd(:)
         real(kind_phys),      intent(in)    :: canopy(:)
         real(kind_phys),      intent(in)    :: tg3(:)
         integer,              intent(in)    :: stype(:)
         real(kind_phys),      intent(in)    :: con_t0c

         integer,              intent(in)    :: nthrds
         character(len=*),     intent(out)   :: errmsg
         integer,              intent(out)   :: errflg

         ! Local variables
         integer :: i, j, ix, vegtyp
         real(kind_phys) :: rsnow

         !--- Noah MP
         integer              :: soiltyp, isnow, is, imn
         real(kind=kind_phys) :: masslai, masssai, snd
         real(kind=kind_phys) :: bexp, ddz, smcmax, smcwlt, dwsat, dksat, psisat

         real(kind=kind_phys), dimension(:), allocatable :: dzsno
         real(kind=kind_phys), dimension(:), allocatable :: dzsnso

         ! Initialize CCPP error handling variables
         errmsg = ''
         errflg = 0

         if (is_initialized) return
         iamin=999
         iamax=-999
         jamin=999
         jamax=-999

!> - Call read_o3data() to read ozone data 
         call read_o3data (ntoz, me, master)

         ! Consistency check that the hardcoded values for levozp and
         ! oz_coeff in GFS_typedefs.F90 match what is set by read_o3data
         ! in GFS_typedefs.F90: allocate (Tbd%ozpl (IM,levozp,oz_coeff))
         if (size(ozpl, dim=2).ne.levozp) then
            write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",    &
                  "levozp from read_o3data does not match value in GFS_typedefs.F90: ", &
                  levozp, " /= ", size(ozpl, dim=2)
            errflg = 1
         end if
         if (size(ozpl, dim=3).ne.oz_coeff) then
            write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",      &
                  "oz_coeff from read_o3data does not match value in GFS_typedefs.F90: ", &
                  oz_coeff, " /= ", size(ozpl, dim=3)
            errflg = 1
         end if

!> - Call read_h2odata() to read stratospheric water vapor data
         call read_h2odata (h2o_phys, me, master)

         ! Consistency check that the hardcoded values for levh2o and
         ! h2o_coeff in GFS_typedefs.F90 match what is set by read_o3data
         ! in GFS_typedefs.F90: allocate (Tbd%h2opl (IM,levh2o,h2o_coeff))
         if (size(h2opl, dim=2).ne.levh2o) then
            write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",     &
                  "levh2o from read_h2odata does not match value in GFS_typedefs.F90: ", &
                  levh2o, " /= ", size(h2opl, dim=2)
            errflg = 1
         end if
         if (size(h2opl, dim=3).ne.h2o_coeff) then
            write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",       &
                  "h2o_coeff from read_h2odata does not match value in GFS_typedefs.F90: ", &
                  h2o_coeff, " /= ", size(h2opl, dim=3)
            errflg = 1
         end if

!> - Call read_aerdata() to read aerosol climatology
         if (iaerclm) then
            ! Consistency check that the value for ntrcaerm set in GFS_typedefs.F90
            ! and used to allocate aer_nm matches the value defined in aerclm_def
            if (size(aer_nm, dim=3).ne.ntrcaerm) then
               write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",     &
                     "ntrcaerm from aerclm_def does not match value in GFS_typedefs.F90: ", &
                     ntrcaerm, " /= ", size(aer_nm, dim=3)
               errflg = 1
            else
               ! Update the value of ntrcaer in aerclm_def with the value defined
               ! in GFS_typedefs.F90 that is used to allocate the Tbd DDT.
               ! If iaerclm is .true., then ntrcaer == ntrcaerm
               ntrcaer = size(aer_nm, dim=3)
               ! Read aerosol climatology
               call read_aerdata (me,master,iflip,idate,errmsg,errflg)
            endif
         else
            ! Update the value of ntrcaer in aerclm_def with the value defined
            ! in GFS_typedefs.F90 that is used to allocate the Tbd DDT.
            ! If iaerclm is .false., then ntrcaer == 1
            ntrcaer = size(aer_nm, dim=3)
         endif

!> - Call read_cidata() to read IN and CCN data
         if (iccn == 1) then
           call read_cidata (me,master)
           ! No consistency check needed for in/ccn data, all values are
           ! hardcoded in module iccn_def.F and GFS_typedefs.F90
         endif

!> - Call tau_amf dats for  ugwp_v1
         if (do_ugwp_v1) then
            call read_tau_amf(me, master, errmsg, errflg)
         endif

!> - Initialize soil vegetation (needed for sncovr calculation further down)
         call set_soilveg(me, isot, ivegsrc, nlunit)

!> - Call setindxoz() to initialize ozone data
         if (ntoz > 0) then
           call setindxoz (im, xlat_d, jindx1_o3, jindx2_o3, ddy_o3)
         endif

!> - Call setindxh2o() to initialize stratospheric water vapor data
         if (h2o_phys) then
           call setindxh2o (im, xlat_d, jindx1_h, jindx2_h, ddy_h)
         endif

!> - Call setindxaer() to initialize aerosols data
         if (iaerclm) then
           call setindxaer (im, xlat_d, jindx1_aer,          &
                            jindx2_aer, ddy_aer, xlon_d,     &
                            iindx1_aer, iindx2_aer, ddx_aer, &
                            me, master)
           iamin = min(minval(iindx1_aer), iamin)
           iamax = max(maxval(iindx2_aer), iamax)
           jamin = min(minval(jindx1_aer), jamin)
           jamax = max(maxval(jindx2_aer), jamax)
         endif

!> - Call setindxci() to initialize IN and CCN data
         if (iccn == 1) then
           call setindxci (im, xlat_d, jindx1_ci,      &
                           jindx2_ci, ddy_ci, xlon_d,  &
                           iindx1_ci, iindx2_ci, ddx_ci)
         endif

!> - Call  cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that can resolve GWs
         if (do_ugwp_v1) then
            call cires_indx_ugwp (im, me, master, xlat_d, jindx1_tau, jindx2_tau,  &
                                  ddy_j1tau, ddy_j2tau)
         endif

         !--- initial calculation of maps local ix -> global i and j
         ix = 0
         do j = 1,ny
           do i = 1,nx
             ix = ix + 1
             jmap(ix) = j
             imap(ix) = i
           enddo
         enddo

         !--- if sncovr does not exist in the restart, need to create it
         if (all(sncovr < zero)) then
           if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters'
           !--- compute sncovr from existing variables
           !--- code taken directly from read_fix.f
           sncovr(:) = zero
           do ix=1,im
             if (landfrac(ix) >= drythresh .or. fice(ix) >= min_seaice) then
               vegtyp = vtype(ix)
               if (vegtyp == 0) vegtyp = 7
               rsnow  = 0.001_kind_phys*weasd(ix)/snupx(vegtyp)
               if (0.001_kind_phys*weasd(ix) < snupx(vegtyp)) then
                 sncovr(ix) = one - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data))
               else
                 sncovr(ix) = one
               endif
             endif
           enddo
         endif

         !--- For RUC LSM: create sncovr_ice from sncovr
         if (lsm == lsm_ruc) then
           if (all(sncovr_ice < zero)) then
             if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM'
             sncovr_ice(:) = sncovr(:)
           endif
         endif
         
         if (errflg/=0) return

         if (iaerclm) then
           call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg)
           if (errflg/=0) return
         end if

         !--- For Noah MP or RUC LSMs: initialize four components of albedo for
         !--- land and ice - not for restart runs
         lsm_init: if (lsm_cold_start) then
           if (lsm == lsm_noahmp .or. lsm == lsm_ruc) then
             if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: initialize albedo for land and ice'
             do ix=1,im
               albdvis_lnd(ix)  = 0.2_kind_phys
               albdnir_lnd(ix)  = 0.2_kind_phys
               albivis_lnd(ix)  = 0.2_kind_phys
               albinir_lnd(ix)  = 0.2_kind_phys
               emiss_lnd(ix)    = 0.95_kind_phys
             enddo
           endif
           if (lsm == lsm_ruc) then
             do ix=1,im
               albdvis_ice(ix)  = 0.6_kind_phys
               albdnir_ice(ix)  = 0.6_kind_phys
               albivis_ice(ix)  = 0.6_kind_phys
               albinir_ice(ix)  = 0.6_kind_phys
               emiss_ice(ix)    = 0.97_kind_phys
             enddo
           endif

           noahmp_init: if (lsm == lsm_noahmp) then
             allocate(dzsno (lsnow_lsm_lbound:lsnow_lsm_ubound))
             allocate(dzsnso(lsnow_lsm_lbound:lsoil)           )
             dzsno(:)    = missing_value
             dzsnso(:)   = missing_value

             tvxy(:)     = missing_value
             tgxy(:)     = missing_value
             tahxy(:)    = missing_value
             canicexy(:) = missing_value
             canliqxy(:) = missing_value
             eahxy(:)    = missing_value
             cmxy(:)     = missing_value
             chxy(:)     = missing_value
             fwetxy(:)   = missing_value
             sneqvoxy(:) = missing_value
             alboldxy(:) = missing_value
             qsnowxy(:)  = missing_value
             wslakexy(:) = missing_value
             taussxy(:)  = missing_value
             waxy(:)     = missing_value
             wtxy(:)     = missing_value
             zwtxy(:)    = missing_value
             xlaixy(:)   = missing_value
             xsaixy(:)   = missing_value

             lfmassxy(:)   = missing_value
             stmassxy(:)   = missing_value
             rtmassxy(:)   = missing_value
             woodxy(:)     = missing_value
             stblcpxy(:)   = missing_value
             fastcpxy(:)   = missing_value
             smcwtdxy(:)   = missing_value
             deeprechxy(:) = missing_value
             rechxy(:)     = missing_value

             snowxy (:)    = missing_value
             snicexy(:,:)  = missing_value
             snliqxy(:,:)  = missing_value
             tsnoxy (:,:)  = missing_value
             smoiseq(:,:)  = missing_value
             zsnsoxy(:,:)  = missing_value
             
             imn          = idate(2)
             
             do ix=1,im
               if (landfrac(ix) >= drythresh) then
                 tvxy(ix)     = tsfcl(ix)
                 tgxy(ix)     = tsfcl(ix)
                 tahxy(ix)    = tsfcl(ix)

                 if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) then
                   tvxy(ix)  = con_t0c
                   tgxy(ix)  = con_t0c
                   tahxy(ix) = con_t0c
                 end if

                 canicexy(ix) = 0.0_kind_phys
                 canliqxy(ix) = canopy(ix)

                 eahxy(ix)    = 2000.0_kind_phys

                 cmxy(ix)     = zero
                 chxy(ix)     = zero
                 fwetxy(ix)   = zero
                 sneqvoxy(ix) = weasd(ix)     ! mm
                 alboldxy(ix) = 0.65_kind_phys
                 qsnowxy(ix)  = zero

!                 if (srflag(ix) > 0.001) qsnowxy(ix) = tprcp(ix)/dtp
                 ! already set to 0.0
                 wslakexy(ix) = zero
                 taussxy(ix)  = zero

                 waxy(ix)     = 4900.0_kind_phys
                 wtxy(ix)     = waxy(ix)
                 zwtxy(ix)    = (25.0_kind_phys + 2.0_kind_phys) - waxy(ix) / 1000.0_kind_phys / 0.2_kind_phys

                 vegtyp       = vtype(ix)
                 if (vegtyp == 0) vegtyp = 7

                 if ((vegtyp == isbarren_table) .or. (vegtyp == isice_table) .or. (vegtyp == isurban_table) .or. (vegtyp == iswater_table)) then

                   xlaixy(ix)   = zero
                   xsaixy(ix)   = zero

                   lfmassxy(ix) = zero
                   stmassxy(ix) = zero
                   rtmassxy(ix) = zero

                   woodxy   (ix) = zero
                   stblcpxy (ix) = zero
                   fastcpxy (ix) = zero

                 else

                   xlaixy(ix)   = max(laim_table(vegtyp, imn),0.05_kind_phys)
!                   xsaixy(ix)   = max(saim_table(vegtyp, imn),0.05)
                   xsaixy(ix)   = max(xlaixy(ix)*0.1_kind_phys,0.05_kind_phys)

                   masslai      = 1000.0_kind_phys / max(sla_table(vegtyp),one)
                   lfmassxy(ix) = xlaixy(ix)*masslai
                   masssai      = 1000.0_kind_phys / 3.0_kind_phys
                   stmassxy(ix) = xsaixy(ix)* masssai

                   rtmassxy(ix) = 500.0_kind_phys

                   woodxy(ix)   = 500.0_kind_phys
                   stblcpxy(ix) = 1000.0_kind_phys
                   fastcpxy(ix) = 1000.0_kind_phys

                 endif  ! non urban ...

                 if (vegtyp == isice_table) then
                   do is = 1,lsoil
                     stc(ix,is) = min(stc(ix,is),min(tg3(ix),263.15_kind_phys))
                     smc(ix,is) = one
                     slc(ix,is) = zero
                   enddo
                 endif

                 snd = snowd(ix)/1000.0_kind_phys  ! go to m from snwdph

                 if (weasd(ix) /= zero .and. snd == zero ) then
                   snd = weasd(ix)/1000.0
                 endif

                 if (vegtyp == 15) then                      ! land ice in MODIS/IGBP
                   if (weasd(ix) < 0.1_kind_phys) then
                     weasd(ix) = 0.1_kind_phys
                     snd       = 0.01_kind_phys
                   endif
                 endif

                 if (snd < 0.025_kind_phys ) then
                   snowxy(ix)   = zero
                   dzsno(-2:0)  = zero
                 elseif (snd >= 0.025_kind_phys .and. snd <= 0.05_kind_phys ) then
                   snowxy(ix)   = -1.0_kind_phys
                   dzsno(0)     = snd
                 elseif (snd > 0.05_kind_phys .and. snd <= 0.10_kind_phys ) then
                   snowxy(ix)   = -2.0_kind_phys
                   dzsno(-1)    = 0.5_kind_phys*snd
                   dzsno(0)     = 0.5_kind_phys*snd
                 elseif (snd > 0.10_kind_phys .and. snd <= 0.25_kind_phys ) then
                   snowxy(ix)   = -2.0_kind_phys
                   dzsno(-1)    = 0.05_kind_phys
                   dzsno(0)     = snd - 0.05_kind_phys
                 elseif (snd > 0.25_kind_phys .and. snd <= 0.45_kind_phys ) then
                   snowxy(ix)   = -3.0_kind_phys
                   dzsno(-2)    = 0.05_kind_phys
                   dzsno(-1)    = 0.5_kind_phys*(snd-0.05_kind_phys)
                   dzsno(0)     = 0.5_kind_phys*(snd-0.05_kind_phys)
                 elseif (snd > 0.45_kind_phys) then
                   snowxy(ix)   = -3.0_kind_phys
                   dzsno(-2)    = 0.05_kind_phys
                   dzsno(-1)    = 0.20_kind_phys
                   dzsno(0)     = snd - 0.05_kind_phys - 0.20_kind_phys
                 else
                   errmsg = 'Error in GFS_phys_time_vary.scm.F90: Problem with the logic assigning snow layers in Noah MP initialization'
                   errflg = 1
                 endif

! Now we have the snowxy field
! snice + snliq + tsno allocation and compute them from what we have

                 tsnoxy(ix,:)  = zero
                 snicexy(ix,:) = zero
                 snliqxy(ix,:) = zero
                 zsnsoxy(ix,:) = zero

                 isnow = nint(snowxy(ix))+1 ! snowxy <=0.0, dzsno >= 0.0

                 do is = isnow,0
                   tsnoxy(ix,is)  = tgxy(ix)
                   snliqxy(ix,is) = zero
                   snicexy(ix,is) = one * dzsno(is) * weasd(ix)/snd
                 enddo
!
!zsnsoxy, all negative ?
!
                 do is = isnow,0
                   dzsnso(is) = -dzsno(is)
                 enddo

                 do is = 1,4
                   dzsnso(is) = -dzs(is)
                 enddo
!
! Assign to zsnsoxy
!
                 zsnsoxy(ix,isnow) = dzsnso(isnow)
                 do is = isnow+1,4
                   zsnsoxy(ix,is) = zsnsoxy(ix,is-1) + dzsnso(is)
                 enddo
!
! smoiseq
! Init water table related quantities here
!
                 soiltyp  = stype(ix)
                 if (soiltyp /= 0) then
                   bexp   = bexp_table(soiltyp)
                   smcmax = smcmax_table(soiltyp)
                   smcwlt = smcwlt_table(soiltyp)
                   dwsat  = dwsat_table(soiltyp)
                   dksat  = dksat_table(soiltyp)
                   psisat = -psisat_table(soiltyp)
                 endif

                 if (vegtyp == isurban_table) then
                   smcmax = 0.45_kind_phys
                   smcwlt = 0.40_kind_phys
                 endif

                 if ((bexp > zero) .and. (smcmax > zero) .and. (-psisat > zero)) then
                   do is = 1, lsoil
                     if ( is == 1 )then
                       ddz = -zs(is+1) * 0.5_kind_phys
                     elseif ( is < lsoil ) then
                       ddz = ( zs(is-1) - zs(is+1) ) * 0.5_kind_phys
                     else
                       ddz = zs(is-1) - zs(is)
                     endif
                     smoiseq(ix,is) = min(max(find_eq_smc(bexp, dwsat, dksat, ddz, smcmax),1.e-4_kind_phys),smcmax*0.99_kind_phys)
                   enddo
                 else                                    ! bexp <= 0.0
                   smoiseq(ix,1:4) = smcmax
                 endif                                   ! end the bexp condition

                 smcwtdxy(ix)   = smcmax
                 deeprechxy(ix) = zero
                 rechxy(ix)     = zero

               endif

             enddo ! ix

             if (errflg/=0) return

             deallocate(dzsno)
             deallocate(dzsnso)

           endif noahmp_init
         endif lsm_init

         is_initialized = .true.

      contains

!
! Use newton-raphson method to find eq soil moisture
!
         function find_eq_smc(bexp, dwsat, dksat, ddz, smcmax) result(smc)
            implicit none
            real(kind=kind_phys), intent(in) :: bexp, dwsat, dksat, ddz, smcmax
            real(kind=kind_phys) :: smc
            real(kind=kind_phys) :: expon, aa, bb, func, dfunc, dx
            integer :: iter
            !
            expon = bexp + 1.
            aa    = dwsat / ddz
            bb    = dksat / smcmax ** expon
            smc = 0.5 * smcmax
            !
            do iter = 1,100
              func  = (smc - smcmax) * aa +  bb * smc ** expon
              dfunc = aa + bb * expon * smc ** bexp
              dx    = func / dfunc
              smc   = smc - dx
              if ( abs (dx) < 1.e-6_kind_phys) return
            enddo
         end function find_eq_smc

      end subroutine GFS_phys_time_vary_init
!! @}

!> \section arg_table_GFS_phys_time_vary_timestep_init Argument Table
!! \htmlinclude GFS_phys_time_vary_timestep_init.html
!!
!>\section gen_GFS_phys_time_vary_timestep_init GFS_phys_time_vary_timestep_init General Algorithm
!! @{
      subroutine GFS_phys_time_vary_timestep_init (                                                 &
            me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, nsswr, fhswr, lsswr, fhour, &
            imfdeepcnv, cal_pre, random_clds,        ntoz, h2o_phys, iaerclm, iccn, clstp,          &
            jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, iflip,            &
            jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm,               &
            jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm,              &
            imap, jmap, prsl, seed0, rann, do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,&
            tau_amf, nthrds, errmsg, errflg)

         implicit none

         ! Interface variables
         integer,              intent(in)    :: me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, &
                                                nsswr, imfdeepcnv, iccn, ntoz, iflip
         integer,              intent(in)    :: idate(:)
         real(kind_phys),      intent(in)    :: fhswr, fhour
         logical,              intent(in)    :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm
         real(kind_phys),      intent(out)   :: clstp
         integer,              intent(in)    :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:)
         real(kind_phys),      intent(in)    :: ddy_o3(:),  ddy_h(:)
         real(kind_phys),      intent(inout) :: ozpl(:,:,:), h2opl(:,:,:)
         integer,              intent(in)    :: jindx1_aer(:), jindx2_aer(:), iindx1_aer(:), iindx2_aer(:)
         real(kind_phys),      intent(in)    :: ddy_aer(:), ddx_aer(:)
         real(kind_phys),      intent(inout) :: aer_nm(:,:,:)
         integer,              intent(in)    :: jindx1_ci(:), jindx2_ci(:), iindx1_ci(:), iindx2_ci(:)
         real(kind_phys),      intent(in)    :: ddy_ci(:), ddx_ci(:)
         real(kind_phys),      intent(inout) :: in_nm(:,:), ccn_nm(:,:)
         integer,              intent(in)    :: imap(:), jmap(:)
         real(kind_phys),      intent(in)    :: prsl(:,:)
         integer,              intent(in)    :: seed0
         real(kind_phys),      intent(inout) :: rann(:,:)

         logical,              intent(in)    :: do_ugwp_v1
         integer,              intent(in)    :: jindx1_tau(:), jindx2_tau(:)
         real(kind_phys),      intent(in)    :: ddy_j1tau(:), ddy_j2tau(:)
         real(kind_phys),      intent(inout) :: tau_amf(:)
         integer,              intent(in)    :: nthrds
         character(len=*),     intent(out)   :: errmsg
         integer,              intent(out)   :: errflg

         ! Local variables
         integer :: i, j, k, iseed, iskip, ix
         real(kind=kind_phys) :: wrk(1)
         real(kind=kind_phys) :: rannie(cny)
         real(kind=kind_phys) :: rndval(cnx*cny*nrcm)

         ! Initialize CCPP error handling variables
         errmsg = ''
         errflg = 0

         ! Check initialization status
         if (.not.is_initialized) then
            write(errmsg,'(*(a))') "Logic error: GFS_phys_time_vary_timestep_init called before GFS_phys_time_vary_init"
            errflg = 1
            return
         end if

         !--- switch for saving convective clouds - cnvc90.f
         !--- aka Ken Campana/Yu-Tai Hou legacy
         if ((mod(kdt,nsswr) == 0) .and. (lsswr)) then
           !--- initialize,accumulate,convert
           clstp = 1100 + min(fhswr/con_hr,fhour,con_99)
         elseif (mod(kdt,nsswr) == 0) then
           !--- accumulate,convert
           clstp = 0100 + min(fhswr/con_hr,fhour,con_99)
         elseif (lsswr) then
           !--- initialize,accumulate
           clstp = 1100
         else
           !--- accumulate
           clstp = 0100
         endif

         !--- random number needed for RAS and old SAS and when cal_pre=.true.
         !    imfdeepcnv < 0 when ras = .true.
         if ( (imfdeepcnv <= 0 .or. cal_pre) .and. random_clds ) then

           iseed = mod(con_100*sqrt(fhour*con_hr),1.0d9) + seed0
           call random_setseed(iseed)
           call random_number(wrk)
           do i = 1,cnx*nrcm
             iseed = iseed + nint(wrk(1)*1000.0) * i
             call random_setseed(iseed)
             call random_number(rannie)
             rndval(1+(i-1)*cny:i*cny) = rannie(1:cny)
           enddo

          do k = 1,nrcm
            iskip = (k-1)*cnx*cny
            do ix=1,im
              j = jmap(ix)
              i = imap(ix)
              rann(ix,k) = rndval(i+isc-1 + (j+jsc-2)*cnx + iskip)
            enddo
          enddo

         endif  ! imfdeepcnv, cal_re, random_clds

!> - Call ozinterpol() to make ozone interpolation
         if (ntoz > 0) then
           call ozinterpol (me, im, idate, fhour, &
                            jindx1_o3, jindx2_o3, &
                            ozpl, ddy_o3)
         endif

!> - Call h2ointerpol() to make stratospheric water vapor data interpolation
         if (h2o_phys) then
           call h2ointerpol (me, im, idate, fhour, &
                             jindx1_h, jindx2_h,   &
                             h2opl, ddy_h)
         endif

!> - Call ciinterpol() to make IN and CCN data interpolation
         if (iccn == 1) then
           call ciinterpol (me, im, idate, fhour,     &
                            jindx1_ci, jindx2_ci,     &
                            ddy_ci, iindx1_ci,        &
                            iindx2_ci, ddx_ci,        &
                            levs, prsl, in_nm, ccn_nm)
         endif

!> - Call  cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that resolve GW-activ
         if (do_ugwp_v1) then
           call tau_amf_interp(me, master, im, idate, fhour, &
                               jindx1_tau, jindx2_tau,       &
                               ddy_j1tau, ddy_j2tau, tau_amf)
         endif
         
!> - Call aerinterpol() to make aerosol interpolation
         if (iaerclm) then
           ! aerinterpol is using threading inside, don't
           ! move into OpenMP parallel section above
           call aerinterpol (me, master, nthrds, im, idate, &
                             fhour, iflip, jindx1_aer, jindx2_aer, &
                             ddy_aer, iindx1_aer,           &
                             iindx2_aer, ddx_aer,           &
                             levs, prsl, aer_nm)
         endif
         
!       Not needed for SCM:
!> - Call gcycle() to repopulate specific time-varying surface properties for AMIP/forecast runs
         ! if (nscyc >  0) then
         !   if (mod(kdt,nscyc) == 1) THEN
         !     call gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit,              &
         !         input_nml_file, lsoil, lsoil_lsm, kice, idate, ialb, isot, ivegsrc,         &
         !         use_ufo, nst_anl, fhcyc, phour, landfrac, lakefrac, min_seaice, min_lakeice,&
         !         frac_grid, smc, slc, stc, smois, sh2o, tslb, tiice, tg3, tref, tsfc,        &
         !         tsfco, tisfc, hice, fice, facsf, facwf, alvsf, alvwf, alnsf, alnwf,         &
         !         zorli, zorll, zorlo, weasd, slope, snoalb, canopy, vfrac, vtype,            &
         !         stype, shdmin, shdmax, snowd, cv, cvb, cvt, oro, oro_uf,                    &
         !         xlat_d, xlon_d, slmsk, imap, jmap)
         !   endif
         ! endif

      end subroutine GFS_phys_time_vary_timestep_init
!! @}

!> \section arg_table_GFS_phys_time_vary_timestep_finalize Argument Table
!! \htmlinclude GFS_phys_time_vary_timestep_finalize.html
!!
!>\section gen_GFS_phys_time_vary_timestep_finalize GFS_phys_time_vary_timestep_finalize General Algorithm
!! @{
      subroutine GFS_phys_time_vary_timestep_finalize (errmsg, errflg)

         implicit none

         ! Interface variables
         character(len=*),                 intent(out)   :: errmsg
         integer,                          intent(out)   :: errflg

         ! Initialize CCPP error handling variables
         errmsg = ''
         errflg = 0

      end subroutine GFS_phys_time_vary_timestep_finalize
!! @}

!> \section arg_table_GFS_phys_time_vary_finalize Argument Table
!! \htmlinclude GFS_phys_time_vary_finalize.html
!!
      subroutine GFS_phys_time_vary_finalize(errmsg, errflg)

         implicit none

         ! Interface variables
         character(len=*),                 intent(out)   :: errmsg
         integer,                          intent(out)   :: errflg

         ! Initialize CCPP error handling variables
         errmsg = ''
         errflg = 0

         if (.not.is_initialized) return

         ! Deallocate ozone arrays
         if (allocated(oz_lat)  ) deallocate(oz_lat)
         if (allocated(oz_pres) ) deallocate(oz_pres)
         if (allocated(oz_time) ) deallocate(oz_time)
         if (allocated(ozplin)  ) deallocate(ozplin)

         ! Deallocate h2o arrays
         if (allocated(h2o_lat) ) deallocate(h2o_lat)
         if (allocated(h2o_pres)) deallocate(h2o_pres)
         if (allocated(h2o_time)) deallocate(h2o_time)
         if (allocated(h2oplin) ) deallocate(h2oplin)

         ! Deallocate aerosol arrays
         if (allocated(aerin)   ) deallocate(aerin)
         if (allocated(aer_pres)) deallocate(aer_pres)

         ! Deallocate IN and CCN arrays
         if (allocated(ciplin)  ) deallocate(ciplin)
         if (allocated(ccnin)   ) deallocate(ccnin)
         if (allocated(ci_pres) ) deallocate(ci_pres)

         ! Deallocate UGWP-input arrays
         if (allocated(ugwp_taulat)) deallocate(ugwp_taulat)
         if (allocated(tau_limb   )) deallocate(tau_limb)
         if (allocated(days_limb  )) deallocate(days_limb)

         is_initialized = .false.

      end subroutine GFS_phys_time_vary_finalize

   end module GFS_phys_time_vary
!> @}
